Mechanisms for Spin-Supersolidity in 5 = 1/2 Spin-Dimer Antiferromagnets 



Mila^ 



00 

o 
o 

(N 



J.-D. Picon, ^'^ A. F. Albuquerque/'"^ K. P. Schmidt,^ N. Laflorencie,'"' M. Troyer,^ and F. 

^ Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland 

''^Institute of Theoretical Physics, Ecole polytechnique federale de Lausanne, Switzerland 

'^School of Physics, The University of New South Wales, Sydney, NSW 2052, Australia 

^ Lehrstuhl fiir theoretische Physik I,Otto-Hahn-StraPe 4, TU Dortmund, D-4-4221 Dortmund, Germany 

"^ Laboratoire de Physique des Solides, Universite Paris-Sud, UMR-8502 CNRS, 91405 Orsay, France 

(Dated: July 25, 2008) 

Using perturbative expansions and the Contractor RenormaUzation (CORE) algorithm, we ob- 
tain effective hard-core bosonic Hamiltonians describing the low-energy physics oi S — 1/2 spin- 
dimer antiferromagnets known to display supersolid phases under an applied magnetic field. The 
resulting effective models are investigated by means of mean-field analysis and Quantum Monte 
Carlo simulations. A "leapfrog mechanism", through means of which extra singlets delocalize in a 
checkerboard-solid environment via correlated hoppings, is unveiled that accounts for the supersolid 
behavior. 
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PACS numbers: 03.75.Nt, 05.30. Jp, 75.10.Jm, 75.40.Mg 

I. INTRODUCTION 

Concepts and techniques developed within a well es- 
tablished research field are often employed in exploring 
new physics displayed by apparently unrelated systems. 
Following this trend, there has been an increased interest 
in field-induced Bose-Einstcin condensation of magnons 
in quantum magnets (for a recent review, see Ref. |lj). Al- 
though the analogy is never complete, this line of research 
undoubtedly has led to considerable success in unveiling 
new phenomena in a growing number of magnetic insu- 
lators under applied magnetic field. The success of this 
approach suggests that one might be able to experimen- 
tally observe more elusive bosonic behavior in quantum 
magnets, such as the phase simultaneously displaying di- 
agonal and off-diagonal order known as supersolid. 

Supersolidity has attracted enormous interest since 
the detection of non-classical rotational inertia in solid 
Helium by Kim and Chanj^!^ Although the correct in- 
terpretation of these measurements is still hotly de- 
bated and there seems to be no consensus on the 
possibility of supersolidity in translationally invariant 
systems^i^i^ii, the occurrence of supersolid phases for 
bosonic models on a lattice is a well established fact. 
While the simplest model of interacting hard-core bosons 
on a square lattice is unstable against phase separa- 
tion, which prevents supersolid behavior f^i^iiS it has been 
shown that frustratio n^^i^^i^^i^"^ , removal of the hard-core 
constrain t "'^^i^^ or inclusion of generalized couplings in the 
Hamiltoniani^iii can stabilize supersolidity. Although a 
more direct implementation of these models remains elu- 
sive, due to the short ranged nature of interactions be- 
tween atoms in optical lattices, one might expect that 
they are relevant in the context of quantum magnets un- 
der an applied magnetic field. 

Indeed, as it was first shown by Ng and Leei^ and fur- 
ther verified by some of us,— an S" = 1/2 spin-dimer 
model on the square lattice with intra-plane coupling 
Ising-likc anisotropy [sec Eq. ([T|) below] has a phase si- 



multaneously displaying diagonal and off-diagonal order, 
the equivalent of a supersolid for spin systems (hence- 
forth dubbed spin-supersolid). Later spin-supersolidity 
was also shown to occur for 5 = 1 systems on a bilayer— 
and on a chain.— However, the exact relationship be- 
tween these spin models and the aforementioned bosonic 
lattice models is not well understood. For instance, one 
might naively expect that the 5 = 1/2 spin-dimer model 
investigated in Refs. [iSl and [l^ will map onto a t — V 
model for hard-core bosons on a square lattice, which is 
known not to display a supersolid phase.— i^ Therefore, 
in order to understand the mechanism behind superso- 
lidity in this model one should analyze the presence of 
extra terms in the effective model. 

Using a perturbative analysis and the contractor renor- 
malization (CORE) method, we derive effective Hamil- 
tonians for the S = 1/2 spin-dimer model studied in 
Refs. [l^ and [l^. A mean-field analysis of the result- 
ing generalized hard-core bosonic Hamiltonian leads to 
a minimal model capable of accounting for supersolid 
behavior, which is then studied by means of Quantum 
Monte Carlo (QMC). 



II. THE MODEL 

We analyze the 5" = 1/2 spin-dimer Hamiltonian ana- 
lyzed by Ng and Lee^^ and some of us^^*^ which reads 



H = J± 2_^ ^i, 



S^. 









i,a = l,2 



(1) 



Si, a is an S" = 1/2 operator attached to the site i of the 
layer a (see Fig. [T]). J± couples spins in different layers 
and is considered to be the essential coupling, being re- 
sponsible for the system's strong dimerized character (we 
set J± = 1 throughout the rest of the paper). Spins in 




FIG. 1: (Color online) (a) The bilayer spin system investi- 
gated in this paper, described by the Hamiltonian of Eq. ([T]). 
The strong coupling Jx, represented by thick vertical lines, 
accounts for the system's strong dimer character. Application 
of a magnetic field along the z direction promotes dimers from 
a singlet (|s), vertical pairs of white circles) to a triplet (|t^), 
pairs of red circles) state and controls the density of emer- 
gent bosons as depicted in (b). Solid (and supersolid) phases 
might be stabilized for field values in the range between the 
lower-critical field /id , where the bottom of the triplet band 
(represented in the inset and separated from the singlet state 
by a zero-field gap 5o and with width D) and the singlet state 
become degenerate, and the upper-critical field /ic2 where the 
system becomes fully polarized. In-plane coupling J leads to 
interactions and hopping amplitudes for emergent bosons, the 
anisotropy A being necessary to stabilize a checkerboard solid 
represented in the upper panel. 



the same layer interact via the couphng J and A is an 
Ising-like anisotropy; finally, the magnetic field h is ap- 
pUed along the easy-axis. We will mainly focus on the set 
of parameters considered in Refs.llSlandlig. J/ J_l = 0.29 
and A = 3.3, leading to an extended supersolid phase as 
evident from the QMC results for the spin stiffness ps 
and static structure factor S'(7r, tt) obtained by Lafloren- 
cie and Mila— and reproduced in Fig. [2l 

Our goal is to show that we can understand the emer- 
gence of SS for the spin model Eq. ([1]) in terms of sim- 
ple microscopic mechanisms. In achieving this, we de- 
rive effective bosonic models for Eq. ([T]) by means of two 
different procedures: high-order perturbative series ex- 
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FIG. 2: (Color online) Quantum Monte Carlo results for the 
spin-dimer model Eq. ([T}. Simulations were performed on a 
16 X 16 X 2 lattice with /3 = 32, for JIJx_ = 1/3.45 « 0.29 
and A = 3.3. (a) Superfluid density ps (open squares) and 
static structure factor 2 x 8(71,77) (filled circles), (b) nor- 
malized magnetization per site m'^ /ml^^ (open diamonds), 
equivalent to particle density in the bosonic language. Dif- 
ferent phases are stabilized as a function of the applied mag- 
netic field, namely: a superfiuid (SF) phase with finite ps 
and vanishing structure factor, an extended supersolid (SS) 
in which both ps and ^(Tr, tt) are finite and a checkerboard 
solid (CBS). Error bars are much smaller than the depicted 
symbols. (Adapted from Ref. [l9l.) 



pansions and the Contractor Renormalization algorithm 
(CORE). We are going to show that correlated hoppings 
for singlets (holes) with amplitudes si (next-nearest- 
neighbor, NNN, hopping which occurs only if at least one 
of the other sites on the same plaquette is occupied) and 
S2 (assisted third-neighbor hopping occurring only when 
the site in between is occupied) , depicted in Fig. [3] (a) and 
(b) respectively, are crucial in accounting for SS behavior 
for the model of Eq. ([1]). It is easy to see [Fig. [3jc)] that 
these processes prevent phase scparatioi*^^ in the hard- 
core bosonic model on the square lattice (t — V model) by 
allowing extra singlets (holes) to delocalize in a checker- 
board solid (CBS) environment by "leapfrogging" on the 
other sublattice and forming a condensate. It is useful to 
define the quantity we call "leapfrog ratio" 
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S2 



1^1 



(2) 



where ti is the nearest-neighbor (NN) hopping amplitude 
for holes. It was shown by Sengupta et alr^ that the 
energetic gain in the domain wall formation behind phase 
separation in the t — V model is cti, where c lies in the 
interval [1,2]. Therefore, for a system of hard-core bosons 
on the square lattice, the energetic gain associated to the 
correlated hoppings depicted in Fig. [3] must be larger 
than cti, implying that the condition S > c/4 must be 
obeyed, for SS behavior to emerge. 
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FIG. 3: (Color online) (a) NNN correlated hopping with am- 
plitude §1 and (b) third-neighbor correlated hopping with am- 
plitude §2' singlets (holes) hop in between red and light-blue 
sites only if the black-filled circles are occupied by holes [in 
(a), at least one of the sites must be occupied; if both are, 
the process occurs with amplitude 2si]. This "leapfrog mech- 
anism" allows for extra holes to delocalize in a checkerboard 
ordered environment, as illustrated in (c), and to condense, 
giving rise to supersolid behavior. 



III. PERTUBATIVE EXPANSIONS 

Following the work of Totsuka^^ and Mila^^ we re- 
strict ourselves to the limit where J± is the main energy 
scale and the system consists of weakly coupled dimcrs. 
The application of a magnetic field lowers the energy of 
one of the triplet bands and at the critical field hd the 
singlet state \s) (holes) and the bottom of the triplet \t^) 
(bosons) band become degenerate [see Fig. (H^b)]. By ex- 
panding the Hamiltonian Eq. ([T]) in terms of the small 
parameter J/J± we can thus obtain an effective hard- 
core bosonic model. Within first-order in J/Jj_, the only 
effective couplings in the model obtained in this way are 
nearest-neighbor hopping amplitude ii and repulsion Vi 
for the emergent bosons (triplets) i^ However, since this 
so-called t — V model [equivalent to Eq. ([3]) below if we set 
0] is known to display no SS phasep^ higher-order 
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effective couplings should be taken into account and we 
proceed to their derivation. 



A. Second-Order Expansion 

We extend the pcrturbative analysis of Mila3^ to sec- 
ond order in J/J±, obtaining the following effective 
Hamiltonian 









{i,j.k) 



(3) 



with effective couplings (we set J^ = 1) 

P2 , J2(2 + A2) ^ 

^^2 = 1 + — ^- — '- + 1 
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FIG. 4: (Color online) A = 3.3. Effective coupUngs ob- 
tained from the PCUTs procedure described in the main text: 
nearest-neighbor hopping amplitude if® (solid dark line) and 
interaction V^^ (dashed dark line) and amplitudes for next- 
nearest-neighbor uncorrelated (^2^, solid light line) and corre- 
lated (s2®, dashed light line) hoppings, as a function of J / J_i_. 
The inset shows the dependence of the ratio sf^/ti"® on J/ J_i_. 



In Eq. ([3]), iii — b\hi = {0, 1} is the occupation number 
for hardcore bosons (|t^) triplets) at the site i on the 
square lattice formed by the spin-dimers. {i,j) denotes 
nearest- neighbor (NN) sites on this lattice and {i,j,k) is 
such that j is a common nearest-neighbor for the sccond- 
or third-neighbor sites i and k^ The physical processes 
at play become more evident after applying a particle- 
hole transformation, (1 — 72;) -^ fii and 6^ _^ 5^ to Eq. ([S]): 



n^i = ^P2 5] n, -f 5] [if (fctfc^. + H.c.) + VT'nSj 



{^,3) 



<! E i^^A + R. 



{i,hk) 



(5) 



^P2. 



We have ignored constant terms and jl = /i -I- 21^/^ 
fii = h\hi is now the singlet (holes) occupation number. 
In addition to the first-order couplings if ^ (NN hopping 
amplitude) and Vi^ (NN repulsion) the second-order ef- 
fective Hamiltonian also contains a correlated hopping 
term with amplitude sf| [see Fig. [Slja-b)].— Correlated 
hoppings have been shown to stabilize supersolidityi^ by 
allowing particles to delocalize in a CBS ordered back- 
ground. However, the second-order amplitude for corre- 
lated hoppings is too smal l^'^i^^ to prevent phase separa- 
tion: the "leapfrog ratio" of Eq. ^, S^^ = 3J/8 « 0.12 
for J/ Jj_ = 0.29 is too small and cannot account for su- 
persolidity. Therefore we extend our analysis and include 
higher-order corrections to the parameters if^, Vj^^ and 
sf I in Eq. ^ with the help of pcrturbative continuous 
unitary transformations (PCUTs). 



B. Perturbative Continuous Unitary 
Transformations (PCUTs) 

The method of continuous unitary transformations 
(CUTs)2^i2^i^i2^ in its perturbative variant22i^i^i2^ and 
quasi-particle conserving form is an efficient tool to de- 
rive effective low-energy models for coupled quantum 
dimer networks in a magnetic field up to high order in 
perturbatioui^i^i 

To this end, the original spin Hamiltonian Eq. ([T]) is 
rewritten in terms of rung triplet operators tj with 
a ~ {±1,0}. This Hamiltonian docs not conserve the 
number of triplets Q = J^i q=±i o ^q^q ^^ ^^^ system. 
The basic idea of quasi-particle conserving CUTs is to 
transform H into Tieff such that [Hcff,Q] — 0, i.e. the 
number of quasi-particles (triplons, in the present case), 
is a conserved quantity. Since the total S^^^ is also a con- 
served quantity, the magnetic field term does not change 
under the unitary transformation. For the case of cou- 
pled dimers in a magnetic field, one can therefore restrict 
to terms in TleS consisting solely of triplet operators t]^ 
in order to describe the low-energy physics. 

In general, a continuous parameter I is introduced such 
that I = refers to the initially given system 7i and 
I = oo corresponds to the final effective system Tieff- Let 
U be the unitary transformation which diagonalizes the 
Hamihonian H and n{l) = UHl)nU{l). Then this uni- 
tary transformation is equivalent to performing an in- 
finite sequence of unitary transforms e^''^'-''^' with the 
anti-hermitian generator 



7?(0 = -u\i)diU{i) 



(6) 



The derivation with respect to / results in the so-called 
flow equation 



diH{l) = [ii{l),7Hl)] 



(7) 



which defines the change of the Hamiltonian during the 
flow. The properties of the effective Haniiltonian depend 
strongly on the choice of the generator i]. Quasi-particle- 
conscrving CUTs chooses -q such that the Hamiltonian 
TYo maps onto an effective Hamiltonian which conserves 
the number of quasi-particlesi^2i2£i2ii22 

In the following we consider the limit of weakly coupled 
rung dimers, i.e. we set J_l = 1 and treat J and A as 
a small expansion parameters. Using a series expansion 
ansatz for 77 and Ti in Eq. [71 one can derive the effective 
quasi-particle conserving Hamiltonian up to high order 
in perturbationi^i^^i^ The results are obtained in the 
thermodynamic limit and in second quantization. 

We stress again that the total S^^^ is a conserved quan- 
tity. The magnetic held term has not changed under 
the unitary transformation and the low-energy physics is 
solely influenced by the local singlet \s) and the triplet 
\t^) polarized parallel to the magnetic fleld (as discussed 
before). Identifying \s) with an empty site and \t^) with 
the presence of a hardcore boson (as before), wc can de- 
duce the effective Hamiltonian in this language by calcu- 
lating matrix elements on finite clustersi^ 



We have extended the derivation of the effective pa- 
rameters appearing in Eq. ([5|), now relabeled as i^®, 
Vi^ and S2^, to sixth order in J / J±^ and have addition- 
ally calculated the amplitude t^^ for uncorrelated third- 
neighbor hopping. Explicit formulas are given in Ap- 
pendix [X] and dependences on J/ J±_ are shown in Fig. [3] 
for A = 3.3. Comparison between t\^ and results ob- 
tained from CORE (see Fig. [S] and discussion below) sug- 
gests that our perturbative analysis remains valid up to 
J/Ji- ^0.15. 

By applying a particle-hole transformation it is pos- 



sible to show that the condition t^^ 



-2s2^, approxi- 



mately fulfilled by our results (Fig. [4]) , implies a vanish- 
ing amplitude for uncorrelated third- neighbor hopping for 
holes (singlets) and therefore the magnitude of s^^ is the 
relevant kinetic scale for supersolidity (a similar situa- 
tion happens for the effective Hamiltonian derived from 
CORE, sec Sec. IV Ap . As we mentioned before, super- 
solid behavior is expected to occur for large enough values 
of S2 ^/t^^ J^iii However, our results for this ratio, shown 
in the inset of Fig. |4l arc clearly too small for preventing 
domain wall formation j-'^Sj — and therefore one does not 
expect to reproduce the extended SS phase observed for 
the original spin model, Eq. ([T]). Consequently, either 
our idea that the model can be described by only taking 
into account \s) and \t^) is wrong, or wc must go beyond 
a perturbative analysis. Since according to Ng and Lee^^ 
contributions from the other two triplets states |i°) and 
|i^^), if non-zero, are negligible close to half-filling, we 
therefore resort on a non-perturbative approach to our 
problem, namely the CORE algorithm. 



IV. CONTRACTOR RENORMALIZATION 

The contractor renormalization (CORE) method was 
introduced by Morningstar and Wcinstein^i^i and has 
been recently^ applied to the study of the spin-dimer 
Hamiltonian described by Eq. ([T|). We extend these re- 
sults by considering the next range in the effective cou- 
plings and analyzing in more detail the resulting effective 
bosonic model. 



A. Procedure 

The basic idea behind CORE (for comprehensive ac- 
counts the reader is referred to Refs. |33 and [401 ) is to 
project out high energy degrees of freedom and to de- 
rive an effective Hamiltonian describing the low-energy 
physics of the original model. Usually this is done by 
first decomposing the lattice on which the original model 
is defined into elementary blocks and diagonalizing the 
Hamiltonian on a single block (while an extended method 
without this restriction was introduced by some of us 
rccentlyf^ here the standard CORE method is more ap- 
propriate). After choosing a suitable number of low- 
energy block states, the model is subsequently diagonal- 
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FIG. 5: Clusters used in the CORE derivation of the effec- 
tive model, labeled according to the longest range effective 
couplings on the square lattice. In this convention, range- 
1 interactions are obtained from the analysis of the cluster 
consisting of two dimers (a), range-2^" from the cluster with 
four dimers forming a square plaquette (b) and range-2 by 
considering three dimers along a line (c). 
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ized on a cluster consisting of a few elementary blocks 
and the lowest energy cluster states are projected onto 
the restricted basis formed by the tensor products of the 
retained block states. An efFective Hamiltonian is then 
obtained by imposing the constraint that the low-energy 
spectrum of the full problem is exactly reproduced and 
by subtracting shorter-range contributions obtained from 
previous steps involving lesser blocks. The validity of 
the procedure can be checked by either analyzing the 
magnitude of long-range efFective couplings (large val- 
ues associated to these signal the inadequacy of the cho- 
sen restricted set of degrees of freedom in accounting for 
the system's low-energy behavior) or, perhaps more ac- 
curately, by keeping track of the weight of the reduced 
density-matrix associated to a single blocki^i^ 

For the spin-dimer model considered here, Eq. ([1]), 
large values for the inter-planc coupling J ^ imply that 
the natural choice when applying CORE is to consider 
the dimers as the elementary blocks: dimer singlet states, 
|s), corresponding to an unoccupied site in the effective 
model living on the square lattice, and an emergent bo- 
son created by promoting one singlet to an S*^ = 1 triplet 
state, \t^)i are the retained block states. The adequacy 
of this reduced set of degrees of freedom in describing 
the Hamiltonian of Eq. ((T|) in the regime know n^^i^^ to 
display supersolid behavior was verified by Abendschein 
and Capponi^ and is confirmed in the present work. 

Our results are obtained from the analysis of the clus- 
ters depicted in Fig. [5l They are labelled according to 
the maximum range for the effective couplings: range-1 
are the results obtained from the analysis of the clus- 
ter containing two dimers shown in Fig. [Slja), range-2^' ^ 
denote the ones from the cluster with four dimers ar- 
ranged as a plaquette [Fig. EJb)] and range-2 results from 
the three-dimcr cluster shown in Fig. [SJc). Wc gauge 
the validity of the mapping onto a system of hard-core 
bosons by analyzing corrections to the nearest-neighbor 
(NN) hopping amplitude ii (for particles) obtained from 
range-2^/^ and range-2 CORE calculations: whenever the 
sum of these contributions exceeds the value obtained 
from range-1 CORE we assume that a valid mapping is 
not obtained. While the criteria used by Abendschein 
and Capponi^ is probably more accurate, our results 



FIG. 6: (Color online) Comparison between CORE (range- 
1, -2^''^ and -2) and PCUTs results for the nearest-neighbor 
hopping amplitude fi (for particles) in the effective bosonic 
model as a function of J/J± for A = 3.3, as in Refs. Ha andll9l. 
The value J/J± — 0.29 is highlighted by the vertical dashed 
line. The vertical solid line indicates the point where longer 
range (2^" and 2) corrections to ti become larger than the 
range-1 contribution, signaling the breakdown of the mapping 
onto a bosonic model (see main text). 



agree qualitatively with theirs and suffice for our analy- 
sis. More importantly, for the parameters (A = 3.3 and 
J/J± = 0.29) leading to supersohdity previously con- 
sidered in the litcraturei^iiS both criteria validate the 
mapping onto the effective bosonic model. 

The efFective hard-core bosonic Hamiltonian obtained 
from the CORE calculation is, after applying a particle- 
hole transformation (1 — n,) -^ fii and 6^ — > 6, given 
by 



ncfi - 



|-/i^ni + 



V, + W, 



% + s^ + n., 



(8) 



where fl'^ is the chemical potential for the holes (singlets) . 
V comprises two-body interactions and W three- and 
four-body interactions; T, S and TZ are the kinetic con- 
tributions: direct and correlated hopping terms. Full ex- 
pressions for each of these terms arc given in Appendix IbI 



B. Comparison with Pertubative Expansion 

Figure \6\ shows our results for the effective nearest- 
neighbor hopping amplitude ii (for particles) obtained 
from CORE (range-1, -2^/^ and -2) and from PCUTs for 
A = 3.3 and as a function of J/J±. As expected, the var- 
ious ranges CORE results agree with the ones obtained 
from PCUTs in the limit of small J/J±, where both re- 
sults are essentially exact. However, for J/J± > 0.15 
higher-order terms in the perturbative expansion start to 
dominate, invalidating the PCUTs analysis. Crucially, 
for the value J/J± = 0.29 considered in Refs. [l^ and 
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TABLE I; Couplings in the effective Hamiltonian obtained 
from CORE [up to range-2, Eqs. lf8l [Bl1lB5|l ] for A = 3.3 and 
J/J± — 0.29. Units are set by J± = 1. 



UM (highlighted by the vertical dashed line Fig. [6]) , the 
PCUTs expansion is clearly invalid, while longer-range 
CORE results are essentially converged. 

These results can be understood if we remark that 
any perturbative expansion about the weakly coupled 
dimer limit is only valid as long as one stays in the zero- 
field rung-singlet phase, with a finite gap to all three 
triplet modes. However, it has been show n^^i^^ that for 
A = 3.3 and J/J± = 0.29 the zero-field ground-state of 
the spin-dimer model Eq. (jT]) displays long-range Neel 
order implying the existence of a quantum critical point 
Jc{h = 0}/J± < 0.29 (evident from poles in Pade analysis 
for the perturbation series) beyond which our perturba- 
tive expansions become meaningless. On the other hand, 
although CORE relies on a strong dimerized character 
(so that dimer singlets and triplets are the relevant local 
degrees of freedom), it does not assume any particular 
ordering and therefore remains valid across the critical 
regime. 
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FIG. 7: (Color online) A = 3.3, J/J± = 0.29. Mean-field 
results for the CBS structure factor S{tt, -k) (filled circles) 
and condensate density po (open squares) as a function of 
the magnetic field h for: (a) the full range-2 CORE effective 
Hamiltonian, Eqs. (|8l IBll - IB5|l . and (b) the minimal model 
of Eq. ((91, with NN hopping amplitude given by Eq. (fTO)l . 
Values for the effective couplings are shown in Table |T] and 
successive phases are labeled as: condensate, supersolid (SS) 
and checkerboard solid (CBS). 



supersolidity in the spin-dimer model of Eq. ([T]) and (b) 
obtaining a simpler effective model amenable to QMC 
simulations (see below) in order to check whether the 
conjectured mechanisms survive after quantum fluctua- 
tions are taken into account. 



MECHANISM FOR SPIN-SUPERSOLIDITY 



Numerical values obtained from CORE for all effec- 
tive couplings (up to rangc-2) appearing in Eqs. ([8l IBIt ^ 
IBSP are shown in Table |l] for the parameters A = 3.3, 
J/J± = 0.29 used in the original QMC simulationsJ^ii^ 
We use the mean-field (MF) approach discussed in Ap- 
pendix [C] and calculate the dependence of the condensate 
density po and CBS order parameter [see Eqs. (|C41IC5|) ] 
on magnetic field h. The results are shown in Fig. [^a). 
The semi-quantitative agreement between these results 
and the QMC data for the original model Eq. ([1]) (shown 
in Fig. [2]) is remarkable if we keep in mind that only con- 
tributions of up to range-2 have been considered in the 
CORE calculation. However, MF approaches are known 
to overestimate supersolid behavio r^'^°i^^ and the effects 
of quantum fluctuations must be carefully analyzed. 

Unfortunately, the effective Hamiltonian obtained 
from CORE, Eqs. ((51 IB1||B5)) . is complex and poses great 
challenges for more unbiased analysis. We therefore use 
the aforementioned MF procedure in gauging the rela- 
tive importance of each term, with a twofold purpose: 
(a) identifying the dominant mechanism accounting for 



Minimal Hamiltonian 



In deciding on a minimal model we should obviously 
take into account the magnitudes associated with each 
term in Eqs. ([8l IBIIIBSP : we start by neglecting all ef- 
fective couplings smaller than O.lif , where if is the NN 
hopping for holes (singlets). Furthermore, since SS takes 
place only close to half-filling, we can also neglect the 
four-body term with coupling W^ [see Eq. (jB2[) ]. The 
resulting model is identical to the second-order effec- 
tive Hamiltonian [Eq. (I5|)],>24 but with strongly renormal- 
ized couplings. In particular, the couplings associated to 
the correlated hoppings sf and S2 [see Fig. [3] (a,b)] are 
considerably larger than predicted by the perturbative 
analysiSf^l as required for SS to emerge. However, the 
MF analysis of the resulting model shows that the ex- 
tra kinetic energy associated to the large effective ampli- 
tudes for correlated hoppings requires the addition of the 
attractive two-body interactions Vf and V^*-' (see Table 
|T| to stabilize a CBS plateau. These considerations lead 



to the minimal model: 



WSin - -M° II "» + H [^T(^-^j + H.c.) + V^n,h, 

+ Y^ |sf [6j(nji +nj2)6fc + H.C. +V^fi,nk^ 
{{i.k}) 



E \^^(3 



{{(hi))) 



h-jbi +H.C. +V^fiihi 



(9) 



Ui = h\hi is the occupation number for holes; {i,j), 
((i,fc)} and {{{i,l)}) denote, respectively, NN, NNN and 
third-NN sites on the square lattice. The correlated hop- 
ping term with amplitude sf [s2 ] is depicted in Fig.[3l^a) 
[Fig. Mh)]- a hole hops between two NNN [third-NN] 
sites i and k [I] only if at least one of their common NN 
sites j'l, j2 [j] is occupied by a hole4^ 

Mean-field results (not shown) for the superfluid den- 
sity ps and the CBS order parameter ^(Tr, n) for the mini- 
mal model of Eq. ([9]), with effective couplings given in Ta- 
ble |I] (for A — 3.3 and J/J± = 0.29), semi-quantitatively 
reproduce the QMC results (shown in Fig. [2]) for the orig- 
inal spin-dimer model, Eq. ([T]). Unfortunately, this pic- 
ture is too simplistic and results from QMC simulations 
(not shown) for this minimal model show that the CBS 
plateau is destroyed by quantum fluctuations, seemingly 
invalidating our analysis. However, the QMC results 
for S{'K, tt) display a rather pronounced peak, indicating 
that our minimal model is close to a borderline where 
the solid phase appears: this is confirmed by the exis- 
tence of an extended CBS plateau (concomitantly with 
a SS phase) in the QMC results obtained by consider- 
ing slightly smaller values for the NN hopping amplitude 
i^ ^ suggesting that terms neglected in the full effective 
model [Eqs. ([8| IBHIBS]) ]. although relatively small, play 
an important role. 

A closer examination of the terms in the full effective 
CORE Hamiltonian [Eqs. (0 iBlllBSl) ] neglected in deriv- 
ing our minimal model Eq. ([9|) shows that the NN corre- 



lated hoppings with amplitudes S3 and Sg [see Eq. (jB4[) 
and Table [I] have exactly the effect of decreasing the 
holes' (singlets') kinetic energy that may stabilize the 
CBS phase. However, the fact that tf and s^, s^ have 
opposite signs also implies that their inclusion in Eq. © 
has the undcsired effect that the resulting minimal model 
would suffer from the sign problem,. In order to circum- 
vent this problem and be able to perform QMC simula- 
tions, we incorporate Sg and Sg in an effective way: we 
notice that in a perfectly ordered CBS background these 
extra hoppings effectively reduce the NN hopping ampli- 
tude if to the value we denote ?"'" given by (to leading 
order) 



ir" = i?-(is 



3 



\S?\ 



(10) 



MF results [Fig.[7]Jb)] for the new minimal effective model 
obtained by the substitution <f -^ i™™ in Eq. © sug- 
gests that the dominant physical processes arc correctly 
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FIG. 8: (Color online) A = 3.3, J/Jj_ = 0.29. Quantum 
Monte Carlo results for the minimal effective model obtained 
from CORE, Eq. (|9]), considering the NN hopping amplitude 
i™'" from Eq. (llOp (values for the couplings are given in Ta- 
ble IT}, for lattice sizes L = 8, 12, 16 and 24. Error bars 
are much smaller than the depicted symbols, (lower panel) 
Superfluid density, ps, (middle panel) CBS order parameter, 
S{Tv,n)/N; (upper panel) singlet density. The temperature 
is set to T = 1/201/ < iT"'/2L (see main text). Successive 
phases are labeled as: superfluid (SF), supersolid (SS) and 
checkerboard solid (CBS). 



taken into account, at least close to half-filling, as we can 
conclude from the excellent agreement with the results 
for the full effective CORE model [Fig. El^a)] 4^ Further- 
more, the SS region visible in Fig. ^h) is expected to 
survive quantum fluctuations, for a sizable leapfrog ratio 
j:{tf") « 0.43 is obtained for A = 3.3 and J/Jj_ = 0.29, 
something confirmed by our QMC simulations below. 



B. Quantum Monte Carlo Simulations 

We have performed QMC simulations, us- 
ing an extended version^ of the ALPS libraries' 
implementation^!^ of the Stochastic Series Expansion 
(SSE) algorithm42i^ We consider the minimal effective 
model of Eq. ([9]) with NN hopping amplitude i™™ given 
by Eq. (|10p 4^ We evaluate the superfluid density ps, 
obtained in terms of the winding numbers Wx and Wy 



1 



PS 



2f3L'- 



■(lOx^ + Wy^) , 



[ID 



where (3 is the inverse temperature and L is the system 
size, and the CBS order parameter 



S{^,^) = jj^{Y,{-ir-'^n,^nr,) 



(12) 



as a function of the magnetic field h. Since we are in- 
terested in accessing ground-state properties, and the 
main kinetic energy scale in the minimal model Eq. ^ is 
/ J± Rd 0.13, we set the temperature to T = l/20i < 



'•1 



Vf™ /2L. It is important to remark that these temper- 
atures are considerably lower than those considered by 
Ng and Lee,— who assumed that J/Jx = 0.29 was the 
relevant energy scale, and this might explain the round 
shape observed in some of their curves. 

QMC results for ps and ^(Tr, tt) for the minimal model 
of Eq. ([ni with NN hopping amplitude given by Eq. ([TU|) . 
using the effective couplings appearing in Table U (A = 
3.3 and J/J±_ = 0.29), are shown in Fig. [i The overall 
agreement with QMC results for the original spin-dimer 
model Eq. ([T])) ^^i^^ shown in Fig. [^ is good and we can 
conclude that the minimal model of Eqs. (|9|10|1 indeed ac- 
curately describes the low-energy physics of the original 
model and, more importantly, that the "leapfrog mecha- 
nism" presented in Sec. [ill is at least partially responsible 
for spin-supersolid behavior. 

Although the just presented results show that the es- 
sential ingredients for spin-supersolidity have been iden- 
tified, it is clear that quantitative agreement is not 
achieved. Specifically, the extent of the SS phase is con- 
siderably smaller in Fig. [5] than in Fig. [51 reversely, the 
CBS phase in the former is about twice as large than 
in the latter. In trying to understand this mismatch it 
is important to keep in mind that supersolidity emerges 
in this model as the result of a delicate balance between 
kinetic and interaction terms. This is evident in the MF 
analysis discussed in Sec. IV Ai which suggests that the 
effective model obtained from CORE is close to a border- 
line and that small variations in the effective couplings 
can have drastic effects. For instance, we have shown 
that the minimal model Eq. ([9|), with effective couplings 
shown in Tablc[ll does not display a CBS phase; however, 
by replacing tf ^ ?]"'" [Eq. dinD] we obtain a CBS phase 
twice as large as expected. 

Therefore, and since the sign-problem precludes us 
from performing QMC simulations for the full effec- 
tive Hamiltonian [Eqs. ([H IBHIB5|l ]. we conjecture that 
that terms ignored in obtaining the minimal model 
[Eqs. (|9|10p ]. even with small couplings, must be included 
in order to better reproduce the results for the original 
model, shown in Fig. [2 Additionally, the NN correlated 
hoppings with amplitudes s^ and s'^ appear to favor SS 
and the fact that we include only their effects in reduc- 
ing tf [Eq. (|10p ] might be responsible for the reduced SS 
phase in Fig. [5J^ Finally, it is not possible to exclude the 
possibility that longer range effective interactions and/or 
neglected triplet excitations {\t^), \t~^)) may be required 
for obtaining quantitative agreement. 



C. Extent of the Supersolid Phase 

We have extended our MF analysis to the full effec- 
tive model Eqs. ([51 IBIHBSP by varying the parameters A 
and J/ J_L, as a function of the magnetic field h. Results 
are shown in Fig. [51 Grey shaded areas represent val- 
ues of A and J/ Jj_ for which no SS phase is stabilized 
within MF for all values of h and only a superfluid and/or 




///, 



FIG. 9: (Color online) Successive phases stabilized for increas- 
ing magnetic field h in the parameter space of the spin-dimer 
model of Eq. ([!]), as obtained from a mean-field analysis of 
the full CORE Hamiltonian, Eqs. P [BT]|B5)) . Grey shaded 
areas correspond to parameters for which no supersolid phase 
is found and in the region marked as "invalid mapping" the 
CORE expansion is invalid (see Sec. IIV.^ . In the remain- 
ing area a supersolid phase is obtained within the mean-field 
approach. The colorful regions indicate parameters for which 
the singlet "leapfrog ratio" [Eq. ([2])] is larger than the thresh- 
old value S(i5"'") > c/4 (with c £ [1,2]), as required for SS 
phases to appear. The cross highlights parameters A = 3.3 
and J/ Ji_ = 0.29 from Refs. llSJ . Il9l . Phases are labeled as: 
Mott insulator (Mo, empty, and Mi, full), superfluid (SF), 
supersolid (SS) and checkerboard solid (CBS). 



CBS phases are obtained. However, SS does appear over 
an extended region in the parameters space within MF. 
Since it is well known that MF tends to overestimate su- 
persolidity, we have also analyzed the "leapfrog ratio" 
E(t""") [see Eq. ([2])] throughout the parameters space. 
As discussed in Sec. [HI the condition S(t™'") > c/4, 
with c S [1,2], must be satisfied for preventing phase 
separation and stabilizing a SS. Regions for which this 
condition is fulfilled are indicated in Fig. [9l we can see 
that the region where the SS phase is likely to occur is 
much smaller than expected from the MF analysis and, 
in particular, no SS is expected within the perturbative 
limit. 



VI. CONCLUSIONS AND OUTLOOK 

Summarizing, we have obtained effective models de- 
scribing the low-energy physics of a spin-dimer model 
[Eq. ([!])] known to exhibit spin-supersolid behavior— li^ 
with the help of perturbative expansions and of the con- 
tractor renormalization (CORE) algorithm. While the 
perturbative analysis, relying on the assumption of a dis- 
ordered ground state with gapped excitations at zero- 
field, does not reproduce the extended supersolid phase 
observed in the original model (Fig[2|), CORE does not 
assume any particular ordering in the system and is 



shown to reproduce the main features obtained from Vi 

more computationally demanding approaches, even when 
a simple mean-field procedure is applied to the obtained 
effective model. 

Furthermore, we identify the mechanism at play be- 
hind spin-supersolidity and we show that the spin- 
supersolid phase exhibited by the S = 1/2 spin-dimer 
model of Eq. ([1]) can be simply understood in terms of 
the "leapfrog mechanism" illustrated in Fig. [31 Basically, 
a sizable amplitude for correlated hoppings allow extra 
holes (singlets) to delocalize on the other sublattice of 
a checkerboard solid, preventing phase separation and ^P6 

leading to supersolid behavior. 

More generally speaking, we are able to describe the 
physics behind complex phenomena in a simple way by 
deriving effective models with only a few terms and 
rather local couplings. The essential physical ingredi- 
ents can be identified even in a low-order perturbative 
analysis, although more sophisticated approaches, such ,/P6 

as PCUTs and CORE, may be required in obtaining the ^ 

effective couplings. We highlight that both PCUTs and 
CORE are immune to the sign problem and can therefore 
be applied to frustrated and fcrmionic systems, some- 
thing which opens interesting research possibilities. 
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The explicit expressions for each term in the effective 
Hamiltonian obtained from CORE, Eq. ([8]), are given 
here (in the expressions below fii = b\hi is the occupa- 
tion number for holes and x, y are unity vectors for the 
square lattice; constant terms arising from applying a 
particle-hole transformation to the bare CORE effective 
Hamiltonian are ignored) . V comprises two-body inter- 
actions 



+ V2 {nifl^ + x + y + flifl^ + x^y) 

+ V^{n.,ni+2x + nihi+2v) , 



(Bl) 



APPENDIX A: PCUTS SIXTH-ORDER 
EFFECTIVE COUPLINGS 



The effective couplings obtained from the PCUTs anal- 
ysis discussed in Sec. IIIIBI are 



and W three- and four-body interactions 

+ W^iriih^+iihi+y + h^+i+y + fii^y + fi^+Si-y)] (B2) 
+ W^ {hifii+sini+yhi+i+y). 



tV' = ^J- 
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The effective single-boson hopping terms in Eq. ([8]) are 



% = tf (6j6,+£ -t- blk+y + B.C.) 
+ i'^{blk+i+y + ifbi+i^y + B.C.) 



(B3) 



10 



Correlated hopping terms are 

Si[blini+x + ni+y)bi+x+y 
+bl{hi+s. + ni^y)bi+x^y + H.c] 
+ S2[blni+xbi+2x + h\ni+yh+2y + H.C.] 

+ s2\b\bi + x{ni + y + fli + x+y + fli-y + fli+x-y) 

+bibi+y{ni+x + fli+x+y + fli^x + fit-x+y) + H.' 

+ S4 [O^bi^xini+yTli^x+y + IT-i—yni+x — y) 

+blbi+y{ni+xni+x+y + ni^xfli^x+y) + H.C.] 

+s^[blbi+x{ni-x + ni+2x) 

+blbi+y{hi^y + n^+2y) + H.C.] 



sC/rt^ 



t^ 



+Sg {b\ni+ini+yhi+x+y + blhi+xni-ybi+x-y + H.c.) , 

(B4) 

and, finally, hoppings simultaneously involving two- 
bosons 



T^i — r^ {b^b-_^^bi+ybi+x+y 



+bibljf^ybi+xbt+x+y + H.c.) 



m. 

+f^{blk+xb,+ybl_^^^^f^ 

+blk+xk^ybl^-_^ + H.C. 



(B5) 



APPENDIX C: MEAN-FIELD PROCEDURE 

Following the Matsubara-Matsuda semiclassical 
approach (^ we write the hard-core boson effective 
models in terms of 5* = 1/2 pseudo-spin variables. We 
start by replacing the commutation relations for bosons 
on the same site i, 



S = (cos (j) sin d, sin (p sin 9, cos 6) which reads 



Hmf — hcff y ,S^ + 2_^ JijSiSj + J.^j {SfSj + SfSj) 



/ J '^i 



(»j> 



E 



Y^ijk'^i \Pj ^k + ^j ^k) + ^ijk^i ^j ^k 



7 , ^ijkl^i ^j Wk^l + ^k^l ) + ^ijkl^i ^j ^kSi 



{i,j,k,l) 



+ 



{i,3,k,l) 



[^'J^ijkli^i Sj S^ Si 



H.C. 



(C3) 



The parameters hcs (one-body) , J (two-body) , K (three- 
body), L (four-body) and M (double exchange) are de- 
fined in terms of the couplings in the effective bosonic 
Hamiltonian. The superscript 2; (-L) accounts for interac- 
tions (hoppings) between sites coupled as in the bosonic 
Hamiltonian. 

In accounting for the different phases of the Hamilto- 
nian Eq. ^ it suffices to consider a site-factorized wave- 
function \tp) ~ Yii l'0i) assuming two-sublatticc long- 
range order (i = A,B). The variational parameters 
{ipA, 4>B, Oa and 9b) are determined by minimizing the 
ground-state energy per site within this subspace. The 
condensate density corresponds in a MF approach to the 
magnetization in the xy plane 



po = -(sin^6'A + sin^6'B) , 

and the CBS structure factor is 

Si-K,!!) = {cos9a- cos9b)^/4: 



(C4) 



(C5) 



In terms of the density of singlets in the sublattice A, 
given by 



riA 



1 + cos 9 a 



(C6) 



[b,,b,]^[blbl]^0 and [&„&!] = 1 (CI) 



by the fermionic anticommutation relations 



{b,,b,} = {bl,bl} = and {b,,bl}^l (C2) 

while retaining the canonical bosonic commutators for 
operators on different sites i, j. This leads to an algebra 
formally equivalent to that of a spin 1/2. 

We then neglect quantum fluctuations by replacing 
the pseudo-spin operators by their mean value, obtain- 
ing a Hamiltonian in terms of classical spins variables 



with a similar definition for the sublattice B (jib), the 
ground state energy per site Eq (up to a constant) for 
the minimal model of Eqs. ((9l fT0)) . 



Eo = 
^V^uaUb 



+ {V^ + V^){nA' + nB') 



(C7) 

(C8) 

4- Ai'^yJnA{l - 7ia)V«s(1 - nB)cos{(j)A - (j^s) (C9) 
+ 2{S^ + 2Sf)nAnB{2-nB-nA) (CIO) 

nA + riB 



+ ih-^L)- 



We can therefore deduce the following trends: 



(Cll) 



11 




0.1 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
"A 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
"A 



FIG. 10; (Color online) Mean-field ground-state energy per site, Eq, for A = 3.3, J/J± — 0.29 and magnetic field h = 2, using 



the effective couplings shown in Table IT] for models comprising the following terms [see ((Qjl] 



(aj ^ 



V^, if, sf and sf , leading 



to a superfluid phase, (b) fi , Vi, Vp , V3 and ff, leading to a CBS phase, (c) The minimal model from Eqs. (f9l llOp . for 
which a spin-supersolid phase is obtained [cf. Fig. [7jb)]. 



v,^ 



and Vf, 



the terms of 



Due to the sign of Vp , 

Eq. dCZl and Eq. (|U8|l favor the CBS because in 
order to minimize them, one must break A — B 
symmetry (cf. Fig. [10] (b)). 



On the contrary, the kinetic term Eq. (jC9p (cf. 
Fig. [10] (a)) favors the SF phase : it is indeed min- 
imal for riA = riB and (j)j^ — (pg ~ tt. There is then 
no symmetry breaking between A and B sublat- 
tices and the latter relation introduces an order in 
the xy plane confirming the presence of a SS phase. 



The last term Eq. (|C10P is more subtle. In the 
case where sf and §2 are negative as presently, 
the contribution of Eq. (|C10|) does not break the 
translational symmetry and therefore only favors 
the SF phase. But it is yet sufficient to induce a SS 
as shown on Fig.[TU] (c). This figure shows that the 



minimization of the total MF energy indeed leads 
to translational symmetry breaking [ua 7^ ns) but 
as the highest density of both is not equal to one, we 
do not obtain a CBS phase but the SS one. What 
is not shown is that minimizing Eq also leads to 
(f'A — 4'B = TT or in other words to an order in the 
xy plane which is the semiclassical equivalent of the 
condensate density. 



For the sake of completeness, let's mention that if we 
would have considered the contributions of s^ and s^ 
indirectly inserted in the minimal Hamiltonian due to 
the sign problem, they would have led to an MF energy 



12 



{nA + nB)cos{<f>A -<I>b) 



In order to minimize it, 



(C12) 
6b must be equal ei- 



ther to or to TT which introduces an order in the xy 
plane. Moreover, ii (J)a — 4>b = t^ (as actually imposed by 
Eq. lC9p . it also leads to a breaking of the A—B symmetry 
as shown by Fig. [TU](c). At the MF level, this correlated 
hopping term alone favors both symmetry breakings con- 
trary to Eq. ((UTOl) . 
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